Violation of particle number conservation in the GW approximation 
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We present a nontrivial model system ol interacting electrons that can be solved analytically in 
the GW approximation. We obtain the particle number from the GW Green's function strictly 
analytically, and prove that there is a genuine violation of particle number conservation if the self- 
energy is calculated non-self-consistently from a zeroth order Green's function, as done in virtually 
all practical implementations. We also show that a simple shift of the self-energy that partially 
restores self-consistency reduces the numerical deviation significantly. 



INTRODUCTION 

Many-body perturbation theory for condensed-matter 
physics allows the Green's function of a system of inter- 
acting electrons to be formulated in a picture of inde- 
pendent quasiparticles moving in an effective potential. 
The key quantity that incorporates the contributions of 
dynamic exchange and correlation to this effective po- 
tential is the self-energy operator S, which in general is 
both nonlocal and energy dependent. It is itself a func- 
tional of the Green's function that can formally be ex- 
pressed through an infinite series of Feynman diagrams. 
In practice, however, approximations for the functional 
form are required, of which the mosij-popular is the GW 
approximation introduced by HedinE This approach re- 
places the infinite series of Feynman diagrams for the 
self-energy operator by a single Fock-like diagram that is 
the product of the Green's function G and the dynami- 
cally screened Coulomb interaction W , hence the name 
of the scheme. Originally the GW approximation was de- 
vised to be applied self-consistently on an equal footing 
with the Hartree and Hartree-Fock approximations, in 
the sense that the Green's function used to generate the 
self-energy be identical to that obtained from it, but the 
computational cost proved prohibitive. Thus GW self- 
energies are traditionally calculated from a zeroth order 
Hartree or LDA Green's function, although a small de- 
gree of self-consistency is sometimes introduced by rigidly 
shifting the zeroth order Green's function on the energy 
axis in such a way that its chemical potential agrees with 
that of the GW Green's function derived from it. It is in 
this fashion that the GW approximation has in the past 
been successfully applied Jjo a wide rangej-pf materials 
including seniiconductorSjaH simple metals^ and transi- 
tion metals The first fully self-consistent GW calcr 
tions for model systems were reported only recentlyjl 
but these consistently showed an undesired loss of struc- 
tural features in|-|the Green's function in disagreement 
with experiment .u 

As future GW calculations for realistic materials are 
therefore unlikely to adopt a fully self-consistent ap- 
proach, the question of particle number conservation, 
long since a subject of debate in this context, has gained 
renewed significance. Since the imaginary part of the 



Green's function is directly related to the density of states 
as well as to the real-space charge density, it must in- 
tegrate to the correct number of electrons contained in 
the system, a condition frequently used to determine the 
numerical accuracy in practical implementations. How- 
ever, while the fully self-consistent GHC,scheme is known 
to conserve the exact particle number,Ll the same, while 
tacitly assumed, has never actually been proven for the 
approximation used in practice, in which the self-energy 
is calculated from a zeroth order Green's function. Nu- 
merical simulations appeared to corroborate this posi- 
tive assumption, with deviations below one percent fully 
within the range expected due to systematic numerical 
errors, and this has already prompted occasional confir- 
mative conjectures in the literature. □ However, in this 
paper we shall show that upon elimination of these er- 
rors, which arise from numerical integration, transfor- 
mation to Fourier space with a finite broadening of the 
quasiparticle peaks, alignment of the chemical potentials 
only within second order perturbation theory, etc., there 
still remains a genuine albeit numerically small viola- 
tion of particle number conservation, as was previously 
demonstrated by accurate molecular calculations for the 
analogous scheme based on time-dependent Hartree-Fock 
rather than Hartree .tiieory that is frequently employed in 
quantum chemistry.ll2l We also demonstrate that restor- 
ing a degree of self-consistency by means of an appropri- 
ate shift in the zeroth order Green's function, which in 
practice is often omitted because of its negligible influ- 
ence on the band structure, further significantly reduces 
the deviation. The significance of this result is that it 
provides a way to minimize the fundamental limit on the 
accuracy with which quantities such as the charge den- 
sity or the total energy can be calculated in the GW 
approximation. 



DESCRIPTION OF MODEL SYSTEM 

For our proof we consider as an analytically solvable 
counter-example a four-site Hubbard cluster with tetra- 
hedral symmetry, populated by two electrons. To the 
best of our knowledge this is in fact the first model system 
forwarded in the literature for which the Green's func- 
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tion and particle number in the GW approximation can 
be calculated strictly analytically, without any numeri- 
cal errors. Its significance for theoretical investigations 
therefore stretches beyond the objective of this paper, 
and for the benefit of the reader we will thus present all 
derivations in sufficient detail. The Hamiltonian is 
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(1) 



where c^^^ and cr(j arc the creation and annihilation 
operators for an electron at site R with spin a, and 
fiRa- = CRg.CR(j. We set e = At and t = 1 while leav- 
ing the onsite interaction strength U variable. To con- 
struct the GW self-energy we start from a zeroth order 
Green's function in the Hartree approximation, noting 
that the exchange-correlation potential in a correspond- 
ing density-functional treatment is a mere constant due 
to spatial symmetries and would not affect the following 
arguments. Spatial symmetries and degeneracy in a also 
require a uniform fractional site occupation of one half 
and accordingly a Hartree potential of U/2 on all sites. 
Analytic diagonalization of the Hartree 4x4 Hamiltonian 
matrix H^j^/ = ei^RR' — t for each spin orientation then 
yields a nondegenerate ground state at zero energy and a 
threefold degenerate excited state at energy e. With the 
electrons in the ground state |0) the zeroth order Green's 
function becomes 



G|r'(^) 



(R|0}(0|RO ^ ^ (R|^)(HRO ^ 
uj — in ^—^ uj — e + iri ' 



where we have omitted the spin index, rj denotes a pos- 
itive infinitesimal. Using the relation (R|0)(0|R') = 1/4 
for the components of the ground-state vector we proceed 
to calculating the polarization propagator in the random- 
phase approximation (RPA), defined through 

^RR' (^) - -2^ /g|r'(^ + ^')G^^'r(^) d^' 



u — e + irj UJ + e — IT] 



(3) 



including a factor 2 for the spin summation. The polar- 
ization propagator is diagonal in eigenvector space, and 
can thus be analytically inverted to yield the screened 
interaction 



R" 
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= f/5RR' + ^E<i^l'^)HR') 



2z 

iy=l 

1 1 



UJ — z + IT] UJ + z — n] 
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with z = [e(e -I- U)]^^"^. A similar analytic expression for 
the screening in a related model system was previously 
given in Ref. The self-energy in the GW approxima- 
tion is a convolution of the zeroth order Green's function 
and the RPA screened interaction, 

^in'H = ^ jG^n'i^ - uj')Wl^Hu^')e^'^' duj' 
3eC/2 (R|0)(0|R') 
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UJ — z — e + irj UJ + z — irj 



(5) 



where S denotes a positive infinitesimal. We have written 
Y^GW such a way as to emphasize that it is diagonal 
in the eigenvectors {v) of the initial Hartree system. On 
the other hand, diagonalization of the full Hamiltonian 
including the self-energy yields the quasiparticle states of 
the interacting electron system, so for our model the two 
are in fact identical, although the corresponding energy 
eigenvalues are not. In this way we can calculate the 
chemical potential fi of the interacting electron system 
exactly from its true quasiparticle properties. By defi- 
nition the chemical potential is identical to the energy 
eigenvalue of the highest occupied quasiparticle state, 
which is zero for the original Hartree system, and at the 
GW level is given implicitly through the self-energy cor- 
rection. 



(6) 



Here we have allowed for the possibility of using a self- 
energy derived from a zeroth order Green's function 
whose chemical potential has been shifted by uj on the 
energy axis. To simulate the effect of a self-consistent 
calculation we determine the shift by requiring the chem- 
ical potentials of the shifted zeroth order Green's func- 
tion and the GW Green's function obtained from it to 
be identical. While this equation is usually solved within 
second order perturbation theory, the simple form of the 
self-energy for our model allows us to derive the exact 
analytic solution 



U 



3eU' 



UJ 
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(7) 



which correctly approaches zero as [/ ^ 0. Accordingly 
we solve Dyson's equation 

G«^(a;) = GH(c.) + GH(c.)E^^(c. - ci)G^^(^) (8) 

in eigenvector space both for uj — and uj as in Eq. (j^), 
and compare the results. In this diagonal form Dyson's 
equation is analytically solvable. For v = Q the self- 
energy contains one pole, which adds a satellite to the 
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quasiparticle peak of the Hartree Green's function and 
yields a quadratic equation for the positions of the poles 
of . Similarly the quasiparticle and satellite struc- 
ture of the V > Q matrix elements is obtained from the 
zeroes of a third order polynomial, reflecting the richer 
spectrum in the self-energy. The GW Green's function 
therefore takes the analytic form 



riGW 



(w -I- Cj) 
(R|0}(0|R') 
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where we have defined the symbols 



(9) 



u 




8z 



(10) 



together with yi = — (^/3) — 2rcos(0/3) and 7/2,3 = 
— (5/3)4- 2rcos[(7r=F(/))/3] as well as the auxiliary quanti- 
ties r = ±\/H, (j) = arccos(g/r3), q = /27 -bc/6 + d/2, 
p — (3c — b'^)/9, and the polynomial coefficients 



-2e + ui- 



U 
4' 



c = -z{z + e) + e[e-^- - \ - -j-, 



a = z(z + e) ( e — — — 
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RESULTS 

From the Green's function in the GW approximation 
(^ the particle number may be obtained by an analytic 
contour integration along a path closed across the upper 
half plane, to sample all occupied states below the chem- 
ical potential. By inspection we note that these are the 
states at xi and j/i. For the total particle number we 
thus obtain 



N 



^ 27ri J 

R 



2{xi~z~e) ^ 6{yi - z - t){yi + z) ^^^^ 

xi - '-- ' 



X2 



{yi -2/2) (2/1 -2/3) 



including a factor 2 for spin summation. If the parti- 
cle number in the GW approximation was conserved, N 
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FIG. 1. Analytically calculated particle number in the GW 
approximation for our two-electron model system as a func- 
tion of the interaction strength U. Applying a rigid shift to 
the zeroth order Green's function such as to align its chem- 
ical potential self-consistently with that of the GW Green's 
function significantly reduces the violation of particle number 
conservation, but fails to remove it completely. 



would have to be a constant with a value of two. In 
particular, it would also have to be independent of the 
interaction strength U, which so far we have not speci- 
fied, but an analysis of the expression ([l^ ) confirms that 
this is not the case. In Fig. |^ we show the calculated par- 
ticle number as a function of J7, both with and without 
applying the shift uj in the zeroth order Green's function. 
In either case the growing deviation from the true value 
as U increases is clearly visible, evidence of a genuine and 
fundamental violation of particle number conservation in 
the GW approximation as usually applied. In absolute 
terms, however, we find that the discrepancy is much re- 
duced when the chemical potentials are aligned in the 
prescribed way, and that it depends only very weakly on 
the interaction strength up to high values of U. For a 
medium correlation of t/ = 4 the numerical deviation 
amounts to an underestimation of merely 0.21% of the 
true particle number. 

A comparison of the weight factors of G"^^ with those 
of the exact Green's function, calculated by numerical 
diagonalization of the Hamiltonian matrix, shows that in 
the absence of the self-consistency shift oj the main error 
in N stems from a serious overestimation of the satellite 
spectrum, while the weight of the quasiparticle peak at 
xi deviates relatively little from the correct value up to 
high correlation strength. When the self-consistency shift 
is applied the satellites are also overestimated, although 
by a lesser amount. In this case, however, this is the 
result of a balanced weight transfer from the quasiparticle 
peaks that has little influence on the integrated spectral 
weight and yields a total particle number in much better 
agreement with the correct value. 
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SUMMARY 

In summary, we have presented a two-electron model 
system for which the particle mimbcr from the Green's 
function in the GW approximation can be derived strictly 
analytically, without any additional inaccuracies that 
have previously beset numerical calculations. Through 
an analysis of this model system we have demonstrated 
that there is a genuine violation of particle number 
conservation in the GW approximation as it is usually 
applied, i.e., with the self-energy calculated non-self- 
consistcntly from a zeroth order Green's function. How- 
ever, we also find that the numerical deviation from the 
exact particle number can be kept low even for an ex- 
tremely strong correlation by introducing a small degree 
of self-consistency in the form of a simple rigid shift of 
the zeroth order Green's function on the energy axis in 
such a way as to align its chemical potential with that 
of the GW Green's function, which is calculated from it 
through Dyson's equation. While the effect of this shift 
on the quasiparticle band structure may be small, our re- 
sults clearly indicate its significance for the calculation of 
quantities such as the particle number and charge density 
that are derived through integration over the complete 
spectral function. 
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